library(apsrtable)
library(plyr)
library(sandwich)
library(apsrtable)
library(lmtest)
library(xtable)

load('survey2.RData')
load('survey3.RData')

###
#News Choice Analysis
###

#Preparing Combined 2019/2020 Study (For News Choice Analysis)
survey3$espntreatment <- NULL
survey3$espn.bias <- NULL
survey3 <- survey3[,names(survey2)]

survey2$espn.headline <- as.character(survey2$espn.headline)
survey3$espn.headline <- as.character(survey3$espn.headline)
survey2$other.headline <- as.character(survey2$other.headline)
survey3$other.headline <- as.character(survey3$other.headline)
survey2$other.source <- as.character(survey2$other.source)
survey3$other.source <- as.character(survey3$other.source)

espn.use.combined <- rbind.data.frame(survey2,survey3)
espn.use.combined$race[which(!(espn.use.combined$race %in% c('Asian or Asian-American','Black or African American','Hispanic or Latino','White')))] <- 'Other'
espn.use.combined$race <- as.character(espn.use.combined$race)
espn.use.combined$age <- as.numeric(as.character(espn.use.combined$age))
espn.use.combined$education <- as.character(espn.use.combined$education)
espn.use.combined$political.interest <- as.character(espn.use.combined$political.interest)
espn.use.combined$political.interest[which(espn.use.combined$political.interest=="")] <- NA
espn.use.combined$gender <- as.character(espn.use.combined$gender)
espn.use.combined$gender[which(espn.use.combined$gender=="")] <- NA

##
#Estimation
##
fit.w.robust <- function(model, cluster.var, dta){	
	robust.se <- function(model, cluster){
		require(sandwich)
		require(lmtest)
		M <- length(unique(cluster))
		N <- length(cluster)
		K <- model$rank
		dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
		uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
		rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
		rcse.se <- coeftest(model, rcse.cov)
		return(list(rcse.cov, rcse.se))
	}
	
	clustervar <- mapply(paste,cluster.var,dta[!(1:dim(dta)[1] %in% model$na.action),cluster.var],sep="")
	vcov <- robust.se(model, clustervar)
	model$se <- sqrt(diag(vcov[[1]]))
	model$coeftest <- vcov[[2]]
	return(model)
}

chose.sports.pooled <- lm(chose.espn ~ pid_3pt + espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined)
chose.sports.pooled <- fit.w.robust(chose.sports.pooled,'rid',espn.use.combined)

chose.sports.sportsonly <- lm(chose.espn ~ pid_3pt + espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined[which(espn.use.combined$coverage.type=='sports'),])
chose.sports.sportsonly <- fit.w.robust(chose.sports.sportsonly,'rid',espn.use.combined[which(espn.use.combined$coverage.type=='sports'),])

chose.sports.politicizedonly <- lm(chose.espn ~ pid_3pt + espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined[which(espn.use.combined$coverage.type=='politicized'),])
chose.sports.politicizedonly <- fit.w.robust(chose.sports.politicizedonly,'rid',espn.use.combined[which(espn.use.combined$coverage.type=='politicized'),])

#Table 4
apsrtable(chose.sports.pooled,chose.sports.sportsonly,chose.sports.politicizedonly,coef.names=c("(Intercept)","Independent","Republican","Perceive ESPN as Liberal","Sports Fan"),model.names=c("All","Sports Only","Politicized Sports Only"),omitcoef=expression(grep(pattern="(Intercept)|pid|espn|sports",coefnames,invert=TRUE)),notes=list("Models condition on demographic variables","Robust standard errors, clustered by Respondent, in parentheses",stars.note),caption="Pr(Select ESPN) by Party, Content Perceptions and Interests")

source.vector <- NA
party.vector <- NA
perception.vector <- NA 
fan.vector <- NA

party.se <- NA
perception.se <- NA 
fan.se <- NA

for(k in 1:length(unique(espn.use.combined$other.source))){
  print(unique(espn.use.combined$other.source)[k])
  chose.sports.pooled <- lm(chose.espn ~ pid_3pt + espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined[which(espn.use.combined$other.source==unique(espn.use.combined$other.source)[k]),])
  chose.sports.pooled <- fit.w.robust(chose.sports.pooled,'rid',espn.use.combined[which(espn.use.combined$other.source==unique(espn.use.combined$other.source)[k]),])
  
  source.vector[k] <- unique(espn.use.combined$other.source)[k]
 
  party.vector[k] <- chose.sports.pooled$coefficients[3]
  perception.vector[k] <- chose.sports.pooled$coefficients[4]
  fan.vector[k] <- chose.sports.pooled$coefficients[5]
  
  party.se[k] <- chose.sports.pooled$coeftest[3,2]
  perception.se[k] <- chose.sports.pooled$coeftest[4,2]
  fan.se[k] <- chose.sports.pooled$coeftest[5,2]
  
  print(summary(chose.sports.pooled))
}

party.frame <- cbind.data.frame(source.vector,party.vector,party.se)
perception.frame <- cbind.data.frame(source.vector,perception.vector,perception.se)
fan.frame <- cbind.data.frame(source.vector,fan.vector,fan.se)

party.frame$lower <- party.frame$party.vector - 2*party.frame$party.se
party.frame$upper <- party.frame$party.vector + 2*party.frame$party.se

perception.frame$lower <- perception.frame$perception.vector - 2*perception.frame$perception.se
perception.frame$upper <- perception.frame$perception.vector + 2*perception.frame$perception.se

fan.frame$lower <- fan.frame$fan.vector - 2*fan.frame$fan.se
fan.frame$upper <- fan.frame$fan.vector + 2*fan.frame$fan.se

#Produce Figure F1
pdf(height=5,width=12,file="figure-f1.pdf")
par(mfrow=c(1,3),mar=c(4.2,10.5,4,1.25))

plot(y=8:1,x=party.frame$party.vector,pch=16,yaxt='n',ylab='',xlab="Coefficient on Republican Partisanship",xlim=c(-.2,.4),main="Partisanship and Pr(Select ESPN)",cex=2)
abline(v=0,lty=2)
axis(side=2,at=8:1,labels=party.frame$source.vector,las=1,cex.axis=1.5)
segments(y0=8:1,x0=party.frame$lower,x1=party.frame$upper,lwd=4)

plot(y=8:1,x=perception.frame$perception.vector,pch=16,yaxt='n',ylab='',xlab="Coefficient on Perceiving ESPN as Liberal",xlim=c(-.2,.4),main="Perceive as Liberal and Pr(Select ESPN)",cex=2)
abline(v=0,lty=2)
axis(side=2,at=8:1,labels=perception.frame$source.vector,las=1,cex.axis=1.5)
segments(y0=8:1,x0=perception.frame$lower,x1=perception.frame$upper,lwd=4)

plot(y=8:1,x=fan.frame$fan.vector,pch=16,yaxt='n',ylab='',xlab="Coefficient on Sports Fandom",xlim=c(-.2,.4),main="Sports Fandom and Pr(Select ESPN)",cex=2)
abline(v=0,lty=2)
axis(side=2,at=8:1,labels=fan.frame$source.vector,las=1,cex.axis=1.5)
segments(y0=8:1,x0=fan.frame$lower,x1=fan.frame$upper,lwd=4)

dev.off()

###
#Marginal Effects of Politicized Attitudes on News Choice by Party
###

marginal.effects <- function(model){
	me.out <- model$coeftest
	point.diff <- me.out[2,1]
	se.diff <- me.out[2,2]
	upper <- point.diff + 2*se.diff
	lower <- point.diff - 2*se.diff
	output <- c(point.diff,se.diff,lower,upper)
	return(output)
}

all.bias <- marginal.effects(fit.w.robust(lm(chose.espn ~ espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined),'rid',espn.use.combined))
dem.bias <- marginal.effects(fit.w.robust(lm(chose.espn ~ espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined[which(espn.use.combined$pid_3pt=='Democrat'),]),'rid',espn.use.combined[which(espn.use.combined$pid_3pt=='Democrat'),]))
ind.bias <- marginal.effects(fit.w.robust(lm(chose.espn ~ espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined[which(espn.use.combined$pid_3pt=='Independent'),]),'rid',espn.use.combined[which(espn.use.combined$pid_3pt=='Independent'),]))
rep.bias <- marginal.effects(fit.w.robust(lm(chose.espn ~ espn.liberal + sports.preference + age + gender + education + race + political.interest,data=espn.use.combined[which(espn.use.combined$pid_3pt=='Republican'),]),'rid',espn.use.combined[which(espn.use.combined$pid_3pt=='Republican'),]))

combined.frame.all <- rbind.data.frame(dem.bias,ind.bias,rep.bias)
names(combined.frame.all) <- c('difference','se','lower','upper')

#Produce Figure 2
pdf(file='figure-2.pdf',height=4,width=6)
par(mar=c(4.1,6.5,1,1.5))
plot(y=3:1,x=combined.frame.all$difference,type='n',ylim=c(.85,3.15),xlim=c(-.15,.15),yaxt='n',ylab='',xlab='Marginal Difference in Choice from Perceiving Liberal Bias at ESPN')
abline(v=0,lty=2)
#points(y=c(4,3,2,1)+.15,x=combined.frame.all$difference,pch=16,cex=2)
#segments(y0=c(4,3,2,1)+.15,x0=combined.frame.all$lower,x1=combined.frame.all$upper,lwd=4)
points(y=c(3,2,1),x=combined.frame.all$difference,pch=16,cex=2,col='black')
segments(y0=c(3,2,1),x0=combined.frame.all$lower,x1=combined.frame.all$upper,lwd=4,col='black')
#points(y=c(4,3,2,1)-.1,x=combined.frame.politicized$difference,pch=16,cex=2,col='gray50')
#segments(y0=c(4,3,2,1)-.1,x0=combined.frame.politicized$lower,x1=combined.frame.politicized$upper,lwd=4,col='gray50')
axis(side=2,at=c(3,2,1),labels=c('Democrats','Independents','Republicans'),las=1)
#legend("topright",pch=c(16,16,16),col=c('black','gray50'),legend=c('Sports','Politicized Sports'),pt.cex=2)
dev.off()
